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Abstract 

p2poc is an add-on toolbox to the Mat lab package pde2path. It is aimed at the numerical 
solution of infinite time horizon optimal control (OC) problems for parabolic systems of PDE 
over ID or 2D spatial domains. The basic idea is to treat the OC problem via the associated 
canonical system in two steps. First we use pde2path to find branches of stationary solutions 
of the canonical system, also called canonical steady states (CSS). In a second step we use the 
results and the spatial discretization of the first step to calculate the objective values of time- 
dependent canonical paths ending at a CSS with the so called saddle point property. This is 
a (typically very high dimensional) boundary value problem (BVP) in time, which we solve by 
combining a modification of the BVP solver TOM with a continuation algorithm in the initial 
states. We explain the design and usage of the package via two example problems, namely the 
optimal management of a distributed shallow lake model, and of a semi-arid grazing system. 
Both show interesting bifurcations of so called patterned CSS, and in particular the latter also a 
variety of patterned optimal steady states. The package (library and demos) can be downloaded 
at WWW. staff.uni-oldenburg.de/hannes.uecker/pde2path. 


1 Introduction 


Denoting the state variable (vector) by n = v{t,x) € and the control by A: = k{t,x) € M, we 
consider spatially distribute infinite time horizon optimal control (OC) problems of the form 


coo 

V{vo{-)):=maxJ{vo{-),k{-,-)), J{vo{-), k{-, ■)) := / e~'^*Jca{vit),k{t))dt, 

Jo 

where Jca{v{-,t),k{-,t)) =-^ [ Jc{v{x,t),k{x,t)) dx 

l“l Jn 


(la) 

(lb) 


is the spatially averaged current value objective function, with the local current value )-M 

a given function; p > 0 is the discount rate, and u : D x [0, oo) —>■ fulfills a PDE of the form 


dtv =-Gi{v,k) := DAv + gi{v,k), v\t=o = vq. (Ic) 

Here, D e is a diffusion matrix, A = + ... + is the Laplacian, and holds in a 

bounded domain f] C with suitable boundary conditions (BC), where for simplicity we restrict 
to homogeneous Neumann djyV = 0, z/ the outer normal. In applications, Jc and Gi of course often 
also depend on a number of parameters, which however for simplicity we do not display herej^ 
Introducing the costates A : 12 x (0, oo) ^ and the (local current value) Hamiltonian 

n = n{v, A, k) = Jc{v, k) + \^{D^v + gi{v, A:)), (2) 

^Gi in 1^) can in fact be of a much more general form, but for simplicity here we stick to &) . The convention 
that dtV — —Gi{v) instead of dtV — Gi{v) with Gi{v) — DAv + gi{v) is inherited via pde2path from the Matlab 
pdetoolbox, which assembles — A into the stiffness matrix K. 
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by Pontryagin’s Maximum Principle (see the references below) for % = ^^7-[(t)dt with the 

spatial integral 

H(t) = / ?{(v(x,t), X(x,t), k(x,t)) dx, (3) 

Jn 

an optimal solution {v^X) (or equivalently (v^k) : Q x [0,oo) ^ has to solve the canonical 

system (CS) 


dfV = dxH = DAv + gi{v,k), = ^o, (4a) 

dtX = pX- dyU = pX + g2{v, k) - DAX, (4b) 

where k = argmax^ T~L{y^ A, fc), which generally we assume to be obtained from solving 

dkH{v^ A, k) = 0. (4c) 

The costates A also fulfill zero flux BC, and derivatives like dyH etc are taken variationally, i.e., 
for T-L. For instance, for A/ = 1 and $(?;,A) := XAv we have $(?;,A) = f^XAvdx = f^(AX)vdx 
by Gaufi’ theorem, hence 6y^{v^ X)[h] = f (AX)hdx, and by the Riesz representation theorem we 
identify A) and hence 9-i;$('L’,A) with the multiplier AA. (Qc) typically applies and yields 

a unique solution k under suitable concavity assumptions on T-L^ and in the absence of control 
constraints, see below. 

In principle we want so solve @ for t G [0, oo), but inj^) we have initial data for only half the 
variables, and in (|^) we have anti-diffusion, such that Qis ill-posed as an initial value problem. 
For convenience we set 

•) •= (^^1^’ .j) : ^ (5) 

and write @ as 

dtu = -G{u, ri) := VAu + f{u), where ^ J)^) > /(“) = 

and where g stands for parameters present, which for instance also includes the discount rate 
p. Besides the boundary condition dj^u = 0 and the initial condition 


v\t =0 = 


(6b) 


we then impose the tranversality condition 


lim e = 0. 

t—)-oo 


(6c) 


A solution u of the canonical system (§ is called a canonical path^ and an equilibrium of ( [6a] ) (which 
automatically fulfils (6c)) is called canonical steady state (CSS). With a slight abuse of notation 
we also call ('u, k) with k given by (§1) a canonical path. See also, e.g., |(^TI15j for more formal 
definitions, further comments on the notions of optimal systems, and, e.g., the significance of the 


transversality condition (6c). 

For general background on OC in a PDF setting see |Trol0j and the references therein, or 
specifically |RZ99a[ IRZ99b| and [A AC 111. Chapters] for Pontryagin’s Maximum Principle for OC 
problems for semilinear diffusive models. However, these works are in a finite time horizon setting, 
and often the objective function is linear in the control and there are control constraints, e.g., 
k{x^t) G K with some bounded interval K. Therefore k is not obtained from the analogue of 
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(&)> but rather takes the values from dK^ which is often called bang control. In, e.g., [CPBT21 
lACKLTl^ . some specific models have been studied in this setting and a rather theoretical way, 
i.e., the focus is on deriving the canonical system and showing well-posedness and the existence of 
an optimal control. |ADS14j additionally contains numerical simulations for a finite time horizon 
control-constrained OC problem for a three species spatial predator-prey system, again leading to 
bang type controls. See also |NPSllj and the references therein for numerical methods for (finite 
time horizon) constrained parabolic optimal control problems. 

Here we do not (yet) consider (active) control or state constraints, and no terminal time, but the 
infinite time horizon. Our models and method are motivated by [BXOSllBXlO] . which also discuss 
Pontryagin’s Maximum Principle in this setting. We do not extend the theory, but rather consider 
after a spatial discretization as a (large) ODE problem, and essentially treat this using the 
notations and ideas from |GCF^08j . to give a numerical framework to calculate optimal solutions. 
Using the canonical system Q we proceed in two steps, which can be seen as a variant of the 
“connecting orbit method”, see, e.g., IBPsnij . and also ^ for further background and remarks 
on the related literature: first we compute (branches of) CSS, and second we compute canonical 
paths connecting to some CSS. This also means that we take a somewhat broader perspective than 
aiming at computing just one optimal control, given an initial condition which without further 
information is ill-posed anyway. Instead, our method aims to give a somewhat global picture by 
identifying the pertinent CSS and their respective domains of attraction. 

(a) CSS branches. We compute (approximate) CSS of (§ , i.e., solutions u of 

G{u,ri)=0, (7) 

together with the BC. For this use the package pde2path |UWB.14l IDR,UW14] to set up a FEM 
discretization of 0 as a continuation/bifurcation problem in one of the parameters, which we call 
7] again. This gives branches r] u{r]) of solutions, which is in particular useful to possibly find 
several solutions j — dX fixed r]. By computing the associated ^^^^) we can 

identify which of these is optimal amongst the CSS. Civen a CSS for simplicity we also write 
Jca{u) := and moreover, have 

J{ii) = Jca{ii)/p> ( 8 ) 

(b) Canonical paths. In a second step (b), we calculate canonical paths ending at a CSS u (and 

often starting at the state values vq of a different CSS uq)^ and the objective values of the canonical 
paths. For this we choose a truncation time T and modify ( [^ to the condition that u{T) G Ws{u) 
and near fi, where Ws{u) denotes the stable manifold of u. In practice, we approximate Ws{u) by 
the stable eigenspace Es{u)^ and thus consider the BVP 

dtu = —G{u)^ (9a) 

v\t =0 = (9b) 

u{T) G Es{u) (and \\u{T) — u\\ small). (9c) 

Using the spatial FEM discretizations, the implementation of G, and the results from the first step 
(a), if the mesh in the FEM consists of n nodes, then u{t) G and ([^) yields a system of 2Nn 

ODEs in the form (with a slight abuse of notation) 

= —G(?i), (10a) 

while the initial and transversality conditions become 


v\t =0 = 

T(?i(T) — ti) = 0 (and \\u{T) — u\\ small). 
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(10b) 

(10c) 




















Here M G jg mass matrix of the FEM mesh, and G (defines the projection 

onto Eu(u). Moreover, @) consists of Nn initial conditions for the states, while the costates A 
(and hence the control k) are free. Thus, to have 2Nn BC altogether we need dimE5(7i) = Nn. On 
the other hand, we always have dim£’ 5 (?i) < Nn^ see |GU15[ Appendix A]. We define the defect 

d{u) := dim£’ 5 ({i) — Nn (11) 


and call a CSS u with d{u) = 0 a CSS with the saddle-point property (SPP). At first sight it may 
appear that d{u) depends on the spatial discretization, i.e., on the number of n of nodes. However, 
d{u) remains constant for finer and finer meshes, see |CU15t Appendix A] for further comments. 

For u = (f), A) with the SPP, and ||'L’o ~ f^ll sufficiently small, we may expect the existence of 
a solution u of which moreover can be found from a Newton loop for ( [Tol ) with initial guess 
u{t) = u. On the other hand, for larger H'L’o ~ a solution of (10) may not exist, or a good initial 
guess may be hard to find, and therefore we use a continuation process also for (10). In the simplest 
setting, assume that for some a G [0,1] we have a solution Ua of (10) with ([Wd) replaced by 


v{0) = avQ + (1 — 


( 12 ) 


ja and use 


as initial guess for 


(e.g., cr = 0 and u = u). We then increase a by some stepsize Sq 
((^), (1^) and ultimately aiming at a = 1. 

To actually solve ([TO^), ([Tol^) and (12) we use TOM |MS02[[MT04[iMSTOQj (see also www.dm. 
uniba. it/~mazzia/mazzia/?page_id=433) in a version mtom which accounts for the mass matrix 
M on the Ihs of (lOr)j^ This predictor {ua) - corrector (mtom for a + 6^) continuation method 
corresponds to the “natural” parametrization of the continuation by a, and is thus implemented 
in p2poc as iscnat (Initial State Continuation NATural). We also give the option to use a secant 
predictor 


U^{t) — ^{t) + 6aT{t)^ r(t) — {u^- 


(j-l)(t) _ ,,,(j-2) 




l(())/||ul 


0 - 1 ) _„ 0 - 2 ) 






(13) 


where and are the two previous steps. However, the corrector still works at fixed a, in 
contrast to the arclength predictor-corrector is care described next. 

It may happen that no solution of ([TO^), ([T^) and (12) is found for a > ao for some ao < 1, 
i.e., that the continuation to the intended initial states fails. In that case, often the BVP shows 
a fold in a, and we use a modified continuation process, letting a be a free parameter and using 
a pseudo-arclength parametrization by a in the BC at t = 0. Since mtom does not allow free 
parameters we add the dummy ODE d = 0, and BCs at continuation step j. 


( 5 , {u{0) — ^^0)) + = cr. 


(14) 


with the solution from the previous continuation step j — 1, and (5,5q;) G x M 

appropriately chosen with ||(s,5q.)||* = 1, where || • ||* is a suitable norm in which may 

contain different weights of v and Va- For s = 0 and Sq. = 1 we find iscnat with stepsize 6^ = cr 
again. To get around folds we may use the secant 

s := (^(?i^-^“^^(0) — ?i^'^“^^(0))/||i^^'^“^^(0) — and Sa — I — ^ 


with small (^, and also a secant predictor 

+ ar (15) 

^mtom is an ad-hoc modification of TOM, and will be replaced by an official version of TOM which handles mass 
matrices once that becomes available. In the following, when discussing, e.g., the behaviour and usage of mtom, we 
note that almost all of this is derived from TOM. 
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for t 1 -^ u^{t) with 

T = ~ u^^~‘^\'))/\\u^^~^\') — u^^~‘^\')\\2 and Tq, = 1 — (^. (16) 

This essentially follows |GCF^08[ §7.2], and is implemented in a routine iscarc (Initial State 
Continuation ARClength). 

Finally, given to calculate T, at startup we solve the generalized adjoint eigenvalue problem 

duG{uf^ = AM$ (17) 

for the eigenvalues A and (adjoint) eigenvectors $, which also gives the defect d{u) by counting the 
negative eigenvalues in A. If d{u) — 0, then from $ E we generate a real base of Eu{u) 

which we sort into the matrix T E 


Acknowledgement. I thank D. Grass, ORCOS Wien, for introducing me to the field of optimal 
control over infinite time horizons, and for clarifying (and posing) many questions regarding the 
aim of software based on Pontryagin’s Maximum Principle for spatially distributed OC problems. 

2 Examples and implementation details 
2.1 The SLOG model 

Following |BX08j . in |GU15j we consider a model for phosphorus P — P{t^x) in a shallow lake, 
and phosphate load k = k{x^t) as a control, which in OD, i.e., in the ODE setting, has analyzed in 
detail for instance in [KWlOj . Here we explain how we set up the spatial so called Shallow Lake 
Optimal Control (SLOG) problem in p 2 poc, and refer to |GU15| for details about the modelling 


and the interpretation of results. The model reads 

roG 

ViPoi-)):=maxJiPo{-),ki;-)), J(Po(-), M', O) := / e-P*JcaiP{t),kit)),dt (18a) 
H-p) Jo 

where Jc{P^ k) = ln(A;) — is the local current value objective function, 

JcaiP{-,t), H-P)) ^ Jc{P{x, t), k{x, t))dx (18b) 

is the spatially averaged current value objective function, and P fulfills the PDE 

dtP{x, t) = k{x, t) - bP{x, t) + + D^Pix, t), (18c) 

1 + P[x,ty 

d^P{x, t)dQ = 0, P{x, t)t=o = Po{x), X E D c (18d) 


The parameter 6 > 0 is the phosphor degradation rate, and 7 > 0 are ecological costs of the 
phosphor contamination P. One wants a low P for ecological reasons, but for economical reasons 
a high phosphate load /c, for instance from fertilizers used by farmers. Thus, the objective function 
consists of the concave increasing function ln(A:), and the concave decreasing function —yP^. We 
consider two scenarios, namely 

Scenario 1: P = 0.5, p = 0.03, y = 0.5,5 E (0.5, 0.8) (primary bif. param.), (19) 

Scenario 2\ D — 0.5, p — 0.3, b = 0.55, y E (2.5, 3.7) (primary bif. param.). (20) 
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With the co-state q and local current value Hamiltonian 


H{P, q, A) = JciP, k) + q[k-bP+ + DAP], 


( 21 ) 


the canonical system for (18) becomes, with k{x^t) = — 


q{x,t)' 


Pix t)'^ 

dtP{x^ t) = /c(x, t) — bP{x^ t) + ^ ^ ^ ^ -h DAP{x^ t), 


1 -h P{x,ty 

dtq{x, t) = 2^P{x, t) + q{x, t) ( p + b- 


2 P(x, t) 


— DAq{x^ t), 


(22a) 

(22b) 

(22c) 

We now explain how to use p2poc to calculate CSS and canonical paths for ( [^ . For this we discuss 
files from the demo directory slocdemo (except for obvious library files), assuming that slocdemo 
is in the same directory as the libraries p2plib, p2poclib and tom. 


(1 + P{x,t)‘^f 

di,P{x,t)ga = 0, d,yq{x,t)ga = 0, P{x,t)t=o ^ Po{x), xEfl. 


2.1.1 Basics of pde2path, and the setup for CSS 

We very briefly review the data structures of pde2path, and refer to |UWR14( IDRUW 1^ for more 
details and the underlying algorithms. The basic structure is a Matlab struct, henceforth called p 
like problem, which has a (large) number of fields (and subfields), as indicated in TableHowever, 


Table 1: Selection (with focus on the semilinear case p. sw. sf em=l) of fields in the structure p describing a 
pde2path problem; see stanparam.m in p2plib for detailed information on the contents of these fields and 
the standard settings. 


field 

fuha 

nc, sw 


purpose 


struct of function handles; in particular the function handles p.fuha.sG, p.fuha.sGjac, 


p.fuha.be, p.fuha.bcjac defining (6a) and Jacobians. 
numerical controls and switches such as p.sw.bifcheck,. 


u,np,nu 

tan,branch 
sol 

usrlam 
eqn,mesh 
plot, file 
mat 


the solution u (including all parameters/auxiliary variables in u(p.nu+l:end)), the number 
of nodes p.np in the mesh, and the number of nodal values p.nu of PDF-variables 
tangent tau(l:p.nu-|-p.nc.nq+l), and the branch, filled via bradat.m and p.fuha.outfu. 
other values/fields calculated at runtime, e.g.: ds (stepsize), res (residual), ... 
vector of user set target values for the primary parameter, default usrlam=[]; 
the tensors c, a, b for the semilinear FEM setup, and the geometry data and mesh, 
switches (and, e.g., figure numbers and directory name) for plotting and file output 
problem matrices, e.g., mass/stiffness matrices M, K for the the semilinear FEM setting. 


most of these can be set to standard values by calling p=stanparam(p). At least in simple problems, 
the user only has to provide: 

1. The geometry of the domain and the boundary conditions. 

2. Function handles (in the semilinear setting of interest here) sG and, for speedup, sGjac, 
implementing G, and its Jacobian. 

3. An initial guess for a solution u of G{u) = 0, i.e., an initial guess for a CSS. 

Typically, the steps 1-3 are put into an init routine, here p=slinit(p,lx,ly,nx,ny,sw,rho), 
where lx,ly,nx,ny are parameters to describe the domain size and discretization, and sw is 
used to set up different initial guesses, see Table The only additions/modifications to the 
standard pde2path setting for CSS problems are as follows: (the additional function handle) 
p.fuha.jc should be set to the local current value objective function, here p.fuha. jc=(§sljcf, 
and p.fuha.outfu to ocbra, i.e., p.fuha.outfu=(§ocbra. This automatically puts Jca{u) at po¬ 
sition 4 of the calculated output-branch. Here we generally use the averaged current objective 
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Table 2: The init routine slinit.m, the rhs slsG.m, the objective function sljcf .m, and the function 
slcon.m. See also, e.g., the source code of slsGjac for the implementation of 


function p=slinit(p,lx,ly,nx,ny,sw) “/o init-routine 

p=stanparam(p) ; “/o set generic parameters to standard, if needed reset below.. 
p.nc.neq=2; p.fuha.sG=@slsG; p.fuha.sGjac=@slsGjac; “/orhs 

p.fuha.outfu=@ocbra; p.fuha.jcf=@sljcf; p.fuha.con=@slcon; % current-val-obj 
p.usrlam=[0.55 0.6 0.65 0.7 0.75]; % target-values for bif-param lam 
[p.mesh.geOjbc]=recnbc2(lx,ly); p.vol=4*lx*ly; % geometry, and volume of dom 
p.fuha.bc=@(p,u) be; p.fuha.bcjac=@(p,u) be; “/o standard Neumann BC 
p.xplot=lx; p. sw. spcalc=0; p.sw.jac=l; p.f ile. smod=100; “/o some more switches 
par=[0.03;0.55;0.5;0.5] ; p.nc.ilam=2; “/o startup param values, and index of main param 
% r=par(l); bp=par(2); cp=par(3); D=par(4); 

p.nc.dsmin=le-6; p.nc.dsmax=0.5; p.nc.lammax=0.8; p.nc.lammin=0.549; p.sol.ds=0.1; 
p=stanmesh(p,nx,ny) ;p=setbmesh(p) ; p. sol.xi=l/p .np; “/o mesh 

p.eqn. c= [1;0;0; 1;-1;0;0;-1] ; p.eqn.a=0; p.eqn.b=0; “/o diffusion tensor and a,b 
switch sw “/o choose initial guess arcoording to switch 

case 1; u=0.3*ones(p.np,1); v=-13*ones(p.np,1); u0=[u v]; p.u=u0(:); % FSC 
case 2; u=2*ones(p.np,1); v=-4*ones(p.np,1); u0=[u v]; p.u=u0(:); % FSM 
case 3; .. % Scenario 2 
end 


p.u=[p.u; par]; p.sw.sfem=l; p=setf emops (p) ; “/o semilin. setting 

[p.u,res]=nloop(p,p.u) ; fprintf ( ^ first res=yog\n^ ,res) ; plotsoKp, 1,1,1) ; 


function r=slsG(p,u) % CS for SLOG, p_t=D*lap p-l/q-b*p+p"2/(l+p"2) 
y q_t=-D lap q+2cp*p+q* (rho+bp-2*p/(l+p''2) "2; 

par=u(p.nu+l:end); r=par(l); bp=par(2); cp=par(3); D=par(4); 

P=u(l:p.np); q=u(p.np+1:2*p.np); 

fl=-l./q-bp*P+P.-2./(l+P.^2); f2=2*cp*P+q.*(r+bp-2*P./(1+P.^2). ^^2); 
f=[fl;f2]; r=D*p.mat.K*u(l:p.nu)-p.mat.M*f; 


function jc=sljcf(p,u) y current value J 

cp=u(p.nu+3: end) ; pv=u(l :p.np) ; kv=-l./u(p .np+1 :p.nu) ; jc=log(kv)-cp*pv.''2; 


function k=slcon(p,u) y extract control from states/costates 
k=-l./u(p.np+l:p.nu); 


function since typically we want to normalize Jd by the domain size for simple comparison between 
different domains. Finally, it is useful to set p. fuha. con=(§slcon, where k=slcon(p,u) extracts 
the control k from the states f, costates A and parameters 77 , all contained in the vector uj^ 

By calling p=cont(p), pde 2 path then first uses a Newton-loop to converge to a (numerical) 
solution, and afterwards attempts to continue in the given parameter. If p. sw.bif check> 0 , then 
pde 2 path detects, localizes and saves to disk bifurcation points on the branch. Afterwards, the 
bifurcating branches can be computed by calling swibra and cont again. These (and other) 
pde 2 path commands (continuation, branch switching, and plotting) are typically put into a script 
file, here bdcmds.m, see Table 

Naturally, there are some modifications to the standard pde 2 path plotting commands, see, e.g., 
plot ID. m. These work as usual by overloading the respective pde 2 path functions by putting the 
adapted file in the current directory. See Fig. [^for example results of running bdemds. 

^Note that we do not use slcon in slsG. However, putting this function for the control into p has the advantage 
that for instance plotting and extracting the value of the control can easily be done by calling some convenience 
functions of p2poc. 
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Figure 1: Example bifurcation diagrams and solution plots from running bdcmds.m. For b < 6foid ~ 0.73 
there are three branches of FCSS, here called FSC (Flat State Clean, low P), FSI (Flat State Intermediate), 
and FSM (Flat State Muddy, high P). On FSCUFSI there are a number of bifurcations to patterned CSS 
branches. The BD in J^a in (b) shows that at, e.g., b = 0.65 we have J(FSI) < J(FMS) < J(pl/pt71) < 
J(pl/ptl6) < J(FSC), suggesting the FSC to be an optimal steady state. However, this does not exclude 
pl/ptl6 or pl/pt71 to be optimal steady states, and moreover, there might be a path from the FSC to some 
other CSS with the SPP, dominating the FSC as a CSS. See |GU15j for further discussion. 

Table 3: Selected commands from the script file bdcmds.m. See the source code for more details, and, e.g., 
bdcmds2D.m for the quite similar commands in 2D. Note that these files are in Matlab cell-modes, and the 
cells should be executed one by one. In the last command, p.fuha.con must be set. 


yo7o stat. BD for sloe, main script-file. First set matlab paths: 
path(^../p2plib^,path);path(^../p2poclib^,path); 
n - Scenario 1, FSC/FSI branch 

close all; p= [];lx=2*pi/0.44; ly=0.1; nx=50; ny=l; sw=l; p=slinit(p,lx,ly,nx,ny,sw); 
p=setfn(p,^f1D; screenlayout(p); p=cont(p,100); 

7o7o -FSM branch 

sw=2; p=slinit(p,lx,ly,nx,ny,sw); p=setfn(p,^f2D; p.nc.dsmax=0.2; p=cont(p,15); 

7o7o -bif from fl (set bpt* and p* and repeat as necessary) 

p=swibra(MlG^bptlG^plG-0.05); p.nc.dsmax=0.3; p.nc.neigdet=50; p=cont(p,150); 

7o7o -plotting of BD, L2 and J_{ca}, and solution plots 

clf(3); pcmp=3; plotbraf ( Ml G ^bptl G3,pcmp, MabG 12, MlG ; 7o FSC (+ other branches) 

clf(3); pcmp=4; plotbraf ( Ml G Mptl G4,pcmp, MabG 12, MlG ; 7o FSC (+other branches) 

plotlDf (^plM ^ptl6Gl,1,1,2) ; plotlDf (^plM ^pt71G2,1,1,2) ; % solution plots ... 
stancssvalf(^plM^ptlGC; 7o extract values <P>, <k>, J_{c,a} from solution in pl/ptl6; 


2.1.2 Canonical paths 

The goal is to calculate canonical paths from some starting state 'L’(O) to a CSS ui with the SPP. 
For this we use one of the continuation algorithm iscnat or iscarc which in turn call mtom, based 
on TOM. Since we only wanted minimal modifications of TOM we found it convenient (though 
somewhat dangerous) to pass a number of parameters to the functions called by mtom via global 
variables. Thus, at the start of the canonical path scripts (here cpdemo.m) we define a number of 
global variables, see Table 

The usage of p2poc to compute canonical paths is best understood by running and inspecting 
the demo file cpdemo.m. Some results of running cpdemo are shown in Fig.[^ (a) shows the “easy” 
case of a canonical path from p3/ptl9 to FSC (up to line 11), while (b) to (f) illustrate the case 
of a fold in a when trying to get a canonical path from pl/pt71 to the FSC. (c)-(e) show the 











































Table 4: Global variables for the computation of canonical paths, i.e., mainly for interfacing the driver 
scripts with the functions called by TOM. 


name 

sO,sl 

Psi 

uO, ul 
par 

uml, um2 
sig 


purpose 

pde2path structs containing the boundary values at t = 0 (sO) and at t = T (si) 

the matrix T to encode the BC at t = T 

vectors containing the current values of u at the boundaries 

the parameter values from si (only for convenience) 

solutions at continuation steps j — 1 and j — 2 (to calculate secant predictors and used in 

extended system in is care) 

current (arclength) stepsize in is care 


two canonical paths obtained from picking two canonical paths from the output of iscarc and a 
posteriori correcting to cr = 0.6 (lines 24-26, only for the “upper” canonical path). Line 21 from 
slcpdemo prepares the calculation of Skiba paths, explained in §2.1.3 We now give a brief overview 
of the involved p 2 poclib functions, with the line-numbers refering to Table 
[alv, vv, sol,udat, tlv,tv,uv] =iscnat (alvin, sol,usee, opt, fn); (line 9) 

Input: alvin as the vector of desired a values, for instance alvin=[0.25 0.5 1 ]. sol,usee 
can be empty (typically on first call), but on subsequent calls should contain the last solution 
and the last secant (if opt.msw=l). fn is a structure containing the filenames for the start 
CSS and end CSS[^ and opt is an options structure containing TOM options and some more, 
see Table [5] and Remark 12.11 

Output: alv as the a vector of successful solves; vv as the canonical path values Ja of the 
successful solves; sol as the (last) canonical path; ydat contains the 2 last steps and the 
last secant (useful for repeated calls, and for using iscnat as startup for iscarc). Finally, 
if opt.retsw=l, then tlv,tv,uv contain data of all successful solves, namely: for the j-th 
step, j = l: length (alv), tlv(j), contains the meshsize (in t), tv(j , 1 : tlv(j)) the mesh, and 
uv (j , 1 : n, 1 : tIv (j) the solution. Thus, via sol. x=xv (j , 1 : tIv (j)); sol. y=squeeze (uv (j , 
l:n, l:tlv(j))); the solution of the j-th step is recovered. This is useful for a posteriori 
inspecting some solution from the continuation (see line 24, and skiba.m in vegdemo). Note 
that uv can be large, and might give memory problems. If opt .retsw= 0 , then tlv,tv,uv are 
empty. 

[alv,vv,usee,esol,tlv,tv,uv] =iscarc(esol,usec,opt,fn); (line 15) 

Input: as for iscnat (without alvin), but with esol containing the extended solution <a), 
and similar for the secant usee; as in iscnat, if opt. start=l, then esol,usee can be empty. 
Output: alv,valG as the vectors of achieved a and usee,esol as the last 

secant/solution; tlv,tv,uv as in iscnat, but all in the sense of extended solutions, i.e. 
(^Uq, , cr). 


Remark 2 . 1 . Concerning the original TOM options we remark that typically we run iscnat and 
iscarc with weak error requirements and what appears to be the fastest monitor and order options, 
i.e., opt .Monitor=3; opt. order=2;. Once continuation is successful (or also if it fails at some cr), 
we can always postprocess by calling mtom again with a higher order, stronger error requirements, 
and different monitor options, e.g., mesh-refinement based on condition rather than error alone. 
See the original TOM documentation. 


The functions iscnat and iscarc are the two main user interface functions for the canonical 

^Taking v from a CSS is basically for convenience: Of course, the initial states vq can be arbitrary, and there 
are no initial conditions for the co-states (and in particular those of the “start CSS” are not used). However, the 
construction is also motivated by the fact that one of the most interesting questions is if given v from some CSS uq 
there exists a canonical path to an end CSS ui, and whether this yields a higher value. Then, it is of course also 
interesting if one can also go the other way round, and for this we provide the flip parameter in setfnflip, see 
below. 
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(a) canonical path from p3/ptl9 to FSC 


(b) A fold in a 




a 


(d) diagnostics 



(e) The “lower” canonical path at a = 0.6 



(f) diagnostics 



Figure 2: Example outputs from running slcpdemo.m 


path numerics. There are a number of additional functions for internal use, and some convenience 
functions, which we briefly review as follows: 

[Psi,mu,d,t] =getPsi(si) ; compute the eigenvalues mu, the defect d, and a suggestion for T. 

Note that this becomes expensive with large 2nN (i.e., the total number of DoF). 

[sol, info] =mtom(0DE,BC,solinit,opt,varargin) ; the ad-hoc modification of TOM, which al¬ 
lows for M in ([Tot)- Extra arguments M and lu,vsw in opt. If opt.lu=0, then \ is used 
for solving linear systems instead of an LU-decomposition, which becomes too slow when 
nN X m becomes too large. See the TOM documentation for all other arguments including 
opt, and note that the modifications in mtom can be identified by searching “HU” in mtom.m. 
Of course mtom (as any other function) can also be called directly (line 26), which for instance 
is useful to postprocess the output of some continuation by changing parameters by hand. 
f=mrhs(t ,u,k); J=fjac(t,u); and f =mrhse (t ,u,k); J=f jace (t ,u) ; the rhs and its Jacobian 
to be called within mtom. These are just wrappers which calculate / and J by calling the 
resp. functions in the pde2path-struct si, which were already set up and used to calculate 
the CSS. si is passed as a global variable. Similar remarks apply to mrhse and f jace for the 
extended setting in iscarc. 

bc=cbcf(ya,yb);[ja,jb]=cbcjac(ya,yb) and bc=cbcfe(ya,yb);[ja,jb]=cbcjace(ya,yb); The 
boundary conditions (in time) for and the associated Jacobians. Implemented by passing 
and similar globally. The *e (as in extended) versions are for iscarc. 
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Table 5: Switches/controls in opt besides the TOM options. Note that iscarc only has a very elementary 
stepsize control via opt. sigmin,opt. sigmax. 


name 

purpose 

name 

purpose 

start 

retsw 

1 for startup, 0 else 

1 to return full continuation data 

M,lu,vsw 

mass matrix, lu-switch and 
(extra) verbosity for mtom 

rhoi 

index of p in par 

tl 

truncation time T 

nti 

7 ^ of points in startup t-mesh 

tv 

current t-mesh 

nsteps 

number of steps for iscarc 

sigmin,sigmax 

min&max stepsize for iscarc 


Table 6: Selected commands from the script file cpdemo.m. See the source code for more details and, e.g., 
more fancy plotting. 


1 “/o driver script for Shallow Lake Optimal Control, first set paths and globals 

2 path( ^ ./tom^path); path( ^ ./p2poclib\path); 

3 close all; clear all; global sO si uO ul Psi par xi yml ym2 sig; 

4 yX Preparations: put filenames into fn, set some bvp parameters 

5 sdO=Ml^ sp0=^ptl2^; sdl=^p3^; spl=^ptl9^ flip=l; I p3->FSC 

6 fn=setfnflip(sdO,spO,sdl,spl,flip); opt=[]; opt=ocstanopt(opt); 

7 opt.rhoi=l; opt.tl=100; opt.start=l; opt.tv=[]; opt.nti=10; opt.retsw=0; 

8 70 % the solve and continue call, and some plots 

9 sol=[]; alvin=[0.1 0.25 0.5 0.75 1]; v=[15,30]; 

10 [alv,vv,sol,ydat,tlv,xv,yv]=iscnat(alvin,sol,[],opt,fn); slsolplot(sol,v); 

11 “/oyo-A fold in alpha, here iscarc needed. Prep, and initial iscarc call 

12 sd0=Ml^; sp0=^ptl2^; sdl=^pl^; spl=^pt71^; flip=l; fn=setfnflip(sd0,sp0,sdl,spl,flip); 

13 esol=[]; ysec=[]; opt.nsteps=3; opt.alvin=[0.2 0.25]; sig=0.1; opt.nti=10; opt.tv=[]; 

14 opt.Stats_step=^on^; opt.start=l; opt.sigmax=l; opt.retsw=l; 

15 [alv,vv,ysecl,esoll,tlv,xv,yv]=iscarc(esol,ysec,opt,fn); opt.start=0; 

16 “/oyo subsequent iccarc-calls (repeat this cell) 

17 opt.nsteps=20; ysec=ysecl; esol=esoll; 7 new input (for repeated calls) 

18 [alvl,vvl,ysecl,esoll,tlvl,xvl,yvl]=iscarc(esol,ysec,opt,fn); 

19 alv=[alv alvl]; vv=[vv vvl]; tlv=[tlv tlvl] ; xv=[xv; xvl] ; yv=[yv; yvl] ; 

20 70% Postprocess sol from iscarc, first a simple plot of J over alpha 

21 alv0=alv; vv0=vv; xv0=xv; yv0=yv; tlv0=tlv; 7 save results for skibademo.m 

22 figure (6); clf; plot (alv(l, :) , vv(l, :) , ^-* O ; xlabeK ^ \alphaO ; ylabel ( M_{a} O ; 

23 7070 fix al from iscarc to some given value and compute CPs, first j=22, then j=34 

24 j=22; tl=tlv(j); n=sl.nu; sol.x=xv(j,1:tlv(j));sol.y=squeeze(yv(j,1:n,1:tlv(j))); 

25 al=0.6; u0=al*s0.u(l:n)+(l-al)*sl.u(l:n); ul=sl.u(l:si.nu); 

26 opt.M=sl.mat.M; sol=mtom(0mrhs,@cbcf,sol,opt); v=[100,30]; slsolplot(sol,v); 


jcaval=jcai(sl,sol,rho) and djca=isjca(sl,sol,rho); 

cT 


J{u) = r 
Jo 


Calculate the objective value 

(23) 

and similarly the normalized 


of the solution u in sol (with Jc taken from sl.fuha.jcf) 
discounted value of a CSS contained in sol.y(: ,end). 
f n=setfnf lip (sd0,sp0,sdl,spl, flip); generate the filename struct fn from sdO, spO (sdO/spO.mat 
contains IC uq) and sdl,spl (contains u); if flip=l, then interchange *0 and *1. 
psol3D(p,sol,wnr,cmp,v,tit); x-t plots of canonical paths; plot component cmp of a canonical 
path sol to figure wnr, with view v and title tit. If cmp=0, then plot the control k, extracted 
from sol via p.fuha.con. 

Thus, after having set up p as in ^2 .1.1 for the CSS, including G and p.fuha. jcf, the user does 
not need to set up any additional functions to calculate canonical paths and their values. However, 
typically there are some functions which should be adapted to the given problem, e.g., for plotting, 
for instance 
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slsolplot(soljv); (line 10) which calls: 

zdia=sldiagn(sol,wnr); Plot some norms on a canonical path as functions of t to figure(wnr). 
This is for instance useful to check the convergence behaviour of the canonical path as t ^ T, 
cf. (d),(f) in Fig.|^ 

2.1.3 A patterned Skiba point 

In ODE OC applications, if there are several locally stable OSS, then often an important issue is to 
identify their domains of attractions. These are separated by so called theshhold or Skiba-points 
(if A = 1) or Skiba-manifolds (if A > 1), see |Ski78| and [GCF^OSi Chapter 5]. Roughly speaking, 
these are initial states from which there are several optimal paths with the same value but leading 
to different CSS. In PDE applications, even under spatial discretization with moderate nA, Skiba 
manifolds should be expected to become very complicated objects. Thus, here we just give one 
example how to compute a patterned Skiba point between ESC and FSM. 

In Line 17-19 of cpdemo.m we attempt to find a path from Pps given by pl/pt71 to (P, g)FSC 
given by FSC/ptl2; this fails due to the fold in a. However, for given a we can also try to find 
a path from the initial state Pa{0) •— <^Pps + (1 ~ <^)^FSC fo the FSM, and compare to the path 
to the ESC. For this, in line 21 of cpdemo.m, we stored the a and Jca values into alvO, vvO, and 
also the path data into tlv0,tv0 uvO. See skibademo.m in Table(in particular line II and the 
following) how to put the values uv0(j , : , 1) into sO and subsequently find the paths to the FSM, 
and Fig. [^for illustration. 

Table 7: The script file skibademo.m. See text for comments. 


1 “/o Skiba example, continues cpdemo.m 

2 % find paths from the yvO initial states from cpdemo.m to FSM 

3 js=10; je=30; jl=js-je+l; % alpha-range; now set the target and Psi to FSM: 

4 sl=loadp(^f2P^ptllD; ul=sl.u(l:si.nu); [Psi,muv,d,tl]=getPsi(sl); 

5 a01=length(alv0); tva=zeros(jl,opt.Nmax+1); % some prep, and fields to hold paths 

6 uva=zeros(jl,n+l,opt .Nmax+1) ; alva= [] ; vva=[]; tavl=[]; sol=[]; 

7 alvin=[0.1 0.25 0.5 0.75 1]; “/o we run from uv0(j,:) to FSM with iscnat 

8 tv=linspace(0,opt.tl,opt.nti); se=2; opt.tv=tv."se./opt.tl''(se-l); doplot=l; 

9 opt.msw=0; opt.Stats_step=^off^; v=[50,8]; % switch off stats 

10 for j=js:je; 

11 fprintf ( ^ j=yoi, al=yog\nP j, alvO(j)); sO .u(l :n)=uv0(j , 1 :n, 1) ^ ; 

12 [alv,vv,soljUdat]=iscnat(alvin,[],[],opt,fn); 

13 if alv(end)==l; Jd=vv0(j)-vv(end) ; fprintf ( Ml-J2=yog\nP Jd) ; % contin. successful 

14 alva=[alva alvO(j)]; vva=[vva vv(end)]; tl=length(sol.x); % put vals in vector 

15 tavl=[tavl tl]; tva(j,1:tl)=sol.x; uva(j,1:n,1:tl)=sol.y; 

16 if abs(Jd)<0.05; doplot=asknu(^plot path?Pdoplot); % Skiba point(s) found 

17 if doplot==l; sol0=[]; alp=alv0(j); % plot the paths to FSC and FSM 

18 solO.x=tv0(j,1:tlvO(j));solO.y=squeeze(uvO(j,1:n,1:tlvO(j))); 

19 psol3Dm(sl,solO,sol, 1,1,[]) ; view(v) ; zlabeK^PD; pause 

20 end 

21 end 

22 end 

23 end 

24 %% plot value diagram 

25 figure(6); plot(alv0(js:jep),vv0(js:jep),^-*bD;hold on;plot(alva,vva,^-*rD; 

26 xlabel(^\alpha^,^FontSizeCsi.plot.fs); ylabel(^J_{a}^,^FontSizeCsi.plot.fs); 
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(a) A Skiba point at a = 0.454 (b) Paths to FSC (blue) and FSM (red) 



Figure 3: Example outputs from skibademo. m 


2.1.4 Further comments 


To keep the demos simple, here we do not include versions of iscnat and is care that use the 
truncation time T as an additional free parameter |GU15^ Fig.5]. However, the files bdcmds.m 
and cpdemo.m also contain some of the commands used to study Scenario 2 , see |GU15|, §2.3] 
for the results, and the directory slocdemo contains the script files bdcmds2D.m and cpdemo2D.m, 
used to compute CSS and canonical paths for ( 22 ) over the 2 D domain D = (—L,L) x (—f, f) 
(based on exactly the same init file slinit.m), some auxiliary plotting functions plotsolf .m and 
plotsolfu.m, and the function sol 2 mov.m used to generate movies of canonical paths. See, e.g., 
|GU15[ Fig. 7,8] for some results, while some movies can be downloaded at the pde 2 path homepage. 


2.2 The vegOC model 


Our second example, from |Uecl5| . concerns the optimal control of a reaction diffusion system used 
to model grazing in a semi arid system for biomass (vegetation) v and soil water u;, following |BX10| . 
Here, semi arid means that there is enough water to support some vegetation, but not enough water 
for a dense homogeneous vegetation. This is an important problem as it is estimated that semi 
arid areas cover about 40% of the world’s land area and support about two billion people, often by 
grazing livestock, www.allcountries.org/maps/world_climate_maps.html, In semi arid areas. 


often overgrazing is a serious threat as it may lead to irreversible desertification, see, e.g., |SBB^09 
and the references therein. 

Denoting the harvesting (grazing) effort as the control by E, we consider 


V{vq,wq) = ma^J{vo,WQ,E), (24a) 

dfV = diAv + [gwp^ — d{l + 6v)]v — Ff, (24b) 

dfW = d2Aw + R{I3 + - {vuV + r^)w, (24c) 


with harvest H = and current value objective function Jc = Jc{v^E) — pH — cE^ which 

thus depends on the price p, the costs c for harvesting/grazing, and E in a classical Cobb-Douglas 
form with elasticity parameter 0 < cr < 1. For the modeling, and the meaning and values of the 
parameters (^, 77 , d, d, /3, g?i, 2 ) we refer to |BX10l[Uecl5| and the references therein (see also 

Table|^ line 9 for the parameter values), and here only remark that p — 0.03,p = 1.1, a — 0.3, c — 
1 are the economic parameters, and we take the rainfall R as the main bifurcation parameter. 
Furthermore, we have the BC and IC 


d^v = diyW = 0 on 5D, {v,w)\t=o = {vq.wq). (24d) 

Denoting the co-states by (A, p) we have the local current value Hamiltonian 

7 /( 7 ;, u;. A, /i, E) = Jc{v^ E) + X [di At; + {gwp^ — d{l + 6v))v — Ff] 

+ p [d2Aw + R{f3 + - {tuV + )«;], (25) 
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and obtain the canonical system 

dtv = T-Lx = di/\v + [gwp^ — d{l + 6v)]v — i7, 
dtw = 7 ^^ = d2Aw + R{^ + ^v) - {tuV + ryj)w, 

dtX = pX — Hv — pX — pav^~^ — X [g{p + l)wv^ — 2dSv — d — av^~^ E^~^] 

— p{R^ — ru)w — diAA, 

dtn = PH - Hw = piJ. - Xgv'^''~^ + p{ruV + r^) - d2Ap, 


where E is obtained from solving Oe'H = 0, giving 


E = 


— XjoL 


V. 


(p-A)(l - a). 

With the notation u = ('L’, le, A,/i), the IC, the BC, and the transversality condition are 
{v^w)\t=o = = 0 on dQ, lim e~^^u{t) = 0. 


t^oo 


(26a) 

(26b) 

(26c) 

(26d) 


(26e) 


(26f) 


To study (26), we write it as dfU = —G{u) and basically need to set up G and the BC. This 


follows the general pde2path settings with the OC related modifications already explained in ^2.1 
and thus we only give the following remarks, first concerning veginit.m, see Table 
In line 2 we only set up p.fuha.sG since in this demo we use p.sw. jac=0 (numerical Jacobians), 
and hence do not need to set p.fuha.sGjac. 

lines 4-7 set the Neumann BC and diffusion tensor for the 4 component system (see gnbc.m and 
isoc.m for documentation) 

lines 8-10 set the desired R values for output of CSS to disk, the parameter values, and the main 
bifurcation parameter. Of course, one could also hard-code all parameters except i?, but 
we generally recommend to treat parameters as parameters since this is needed if later a 
continuation in some other parameter is desired, and since it usually makes the code more 
readable. 


Table 8: Selected commands from the init-routine veginit.m. See the source code for more details. 


1 function p=veginit(p,lx,ly,nx,ny,sw,rho) % init-routine for vegOC 

2 p=stanparam(p); p.nc.neq=4; p.fuha.sG=@vegsG; p.fuha.jcf=@vegjcf; p.fuha.outfu=@ocbra; 

3 p.mesh.geo=rec(lx,ly); p=stanmesh(p,nx,ny); p.sol.xi=0.005/p.np; /o generate mesh 

4 q=zeros(p.nc.neq); g=zeros(p.nc.neq,1); /o setting up Neumann BC for 4 components 

5 bc=gnbc(p.nc.neq,4,q,g); p.fuha.bc=@(p,u) be; p.fuha.bejac=@(p,u) be; 

6 p.dl=0.05; p.d2=10; p.eqn.a=0; p.eqn.b=0; /o setting up K for 4 components 

7 c=diag([p.dl, p.d2, -p.dl, -p.d2]); p.eqn.c=isoc(c,4,1); 

8 p.usrlam=[4 10 20 26 28] ; /o desired R values for output of CSS 

9 par=[rho le-3 0.5 0.03 0.005 0.9 le-3 34 0.01 0.1 1 1.1 0.3]; % par-values 

10 p.nc.ilam=8; % choose the active par, here Rainfall R 

11 /o now continue with setting a few more param and the initial guess ... see veginit.m 


Table 1^ shows the complete codes for setting up G and the current value Jc> Both use the auxiliary 
function efu, which is also used in, e.g., valf .m to tabulate characteristical values of CSS. 

Figure shows a basic bifurcation diagram of CSS in ID with D = (—L,L), L = 5, from the 
script file vegbdld.m, which follows the same principles as the one for the SLOC demo. The blue 
branch in (a) represents the primary bifurcation of PCSS, which for certain R have the SPP, and, 
moreover, are POSS. See also |Uecl5j for more plots, including a comparison with the uncontrolled 
case of so called “private optimization”, and 2D results for D = (—L,L) x (—v^L/2, v^L/2) 
yielding various POSS, including hexagonal patterns. 
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Table 9: Implementation of the rhs G and Jc for (26), and the aux function efu. 


1 function r=vegsG(p,u) % rhs for vegOC problem 

2 par=u(p.nu+l: end) ; rho=par(l) ; g=par(2) ; eta=par(3) ; “/o extract param 

3 d=par(4); del=par(5); beta=par(6); xi=par(7); rp=par(8); up=par(9); rw=par(10); 

4 cp=par(ll); pp=par(12); al=par(13); [e,h, J] =efu(p,u); “/o calculate H 

5 v=u(l:p.np); w=u(p.np+1:2*p.np); % extract soln-components, states 

6 ll=u(2*p.np+1:3*p.np); 12=u(3*p.np+l:4*p.np); % co-states 

7 fl=(g*w.*v."eta-d*(l+del*v)).*v-h; f2=rp*(beta+xi*v)-(up*v+rw).*w; % fl,f2 

8 f3=rho*ll-pp*al*h./v-11.*(g*(eta+1)*w.*v."eta-2*d*del*v-d-al*h./v)-12.*(rp*xi-up*w); 

9 f4=rho*12-ll .* (g*v(eta+1) )-12 .* (-up*v-rw) ; f = [f 1; f 2; f 3; f 4] ; 

10 r=p.mat.K*u(l:p.nu)-p.mat.M*f; % the residual 

1 function jc=vegjcf(p,u); [e,h,jc]=efu(p,u); % J_c for vegOC, here just an interface 

1 function [e,h,J] =efu(p,varargin) “/o extract [e,h,J] from p or u 

2 if nargin>l u=varargin{l}; else u=p.u; end 

3 par=p.u(p.nu+1:end); cp=par(ll); pp=par(12); al=par(13); 

4 v=u(l:p.np); ll=u(2*p.np+l:3*p.np); 

5 gas=( (pp-11) * (1-al) ./cp) . " (1/al) ; e=gas.*v; h=v.''al. *e. " (1-al); 

6 J=pp*v. "al. *e.(1-al)-cp*e; 


(a) 


(b) 


(c) p1/pt11,J_,^=22 




Figure 4: Outputs of bdlddemo .m. (a),(b) bifurcation diagrams of CSS in ID; (c) example solutions. 

The script files vegcpdemo.m for canonical paths, and vegskiba.m for a Skiba point between 
the flat optimal steady state FSS/ptl3 and the POSS pl/pt34, again follow the same principles 
as in the the SLOC demo. See Figures]^ and for example outputs, and |Uecl5j for a detailed 
discussion. In a nutshell, we find that: 

(a) For large R the FCSS is the unique CSS of ( [^ , and is optimal, hence a globally stable FOSS 
(Flat Optimal Steady State). 

(b) For smaller R there are branches of (locally stable) POSS (Patterned Optimal Steady States), 
which moreover dominate all other CSS. 

(c) For the uncontrolled problem, Flat Steady States (FSS) only exist for much larger R than 
the FCSS under control. 

(d) At equal i?, the profit J (or equivalently the discounted value Jdp) of the uncontrolled FSS 
is much lower than the value of the FCSS under control. 
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(a) R = 26, the canonical path from the lower PCSS(pl/ptll) to the FCSS (FSS/ptl7) 


p1/pt11 to FSS/pt17 
Jq=733.37, J=735.41, J^=748.29 



t 



(b) R = 26, the canonical path from the FCSS to the upper PCSS (pl/pt38) 



Figure 5: Example output of vegcpdemo.m: Two canonical paths. The leftmost panels indicate the conver¬ 
gence behaviour, the current value profits, and obtained objective values. The middle and right panels show 
the strategy E and the corresponding behaviour of (v^w). See |Uecl5j for comments and more details. 


(a) A Skiba point at a 


0.9 (b) Paths of (almost) equal values to the FCSS and the upper PCSS. 



Figure 6: R = 28, example outputs from skibademo .m. In (a), the blue line gives J for the canonical path 
t u{t) from (i;,re)Q,(0) := a{v,w)ps + (1 ~ '^)fsS 5 where FSS denotes FSS/ptl3, and PS denotes 

pl/ptl6. The red line gives J for the CP t i-^ u{t) from Pa{0) to the upper PCSS pl/pt34. Similarly, the 
white surfaces in (b) are for u and the colored ones for ft. R = 2S. 


3 Summary and outlook 


With p2poc we provide a toolbox to study OC problems of class Q in a simple and convenient 
way, in ID and 2D. The class 0 is quite general, and with the pde2path machinery we have 
a rather powerful tool to study the bifurcations of CSS. The computation of canonical paths is 
comparatively more involved. Essentially, our step (b) implements for the class 0 (parts of) 
the methods explained for ODE problems in GCF^O^ Chapter 7], and implemented in OCMat 
orcos.tuwien.ac.at/research/ocmat software/, see also [Gra m for an extension of OCMat to 
ID systems of class 0 . In a somewhat more general sense, step (b) is a special case (for PDEs) of 
the “connecting orbit method”. See |DCF^97[ IBPS01| and the references therein for earlier work 
on connecting orbits in ODE problems, including connecting orbits to periodic solutions, which 
for ODE OC problems may also be important as long-run optimal solutions, again cf. [GCF^08| . 
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Our setup for (b) is reasonably fast for up to 4000 degrees of freedom of u at fixed time, e.g., 1000 
spatial discretization points and 4 components, and up to 200 temporal discretization points, i.e., 
up to these values a continuation step in the calculation of a canonical path takes up to a few 
minutes on a desktop computer. 

Of course, there is a rather large number of issues we do not address (yet). Besides periodic 
long-run optimal solutions, one of these are state or control inequality constraints that frequently 
occur in OC problems. For instance, in the SLOG model we need non-negativity of P and fc, and 
similarly oi v^w and E in the vegOC model. In our examples we simply checked these a posteriori 
and found them to be always fulfilled, i.e., inactive. If such constraints become active the problem 
becomes much more complicated. Some extensions in this direction will be added as required by 
examples. 

Clearly, it is tempting to recombine steps (a) and (b) again, at least for specific purposes. One 
example would be the continuation of canonical paths in a parameter rj. Naively, this could be 
done “by hand” by using a canonical path between vq (or vo{r])) and u{r]) as an initial guess 

for a canonical path u{-^r]-\- 6) between vq (or vo{r] -\- 5)) and u{r] + 5), at the parameter value r]-\-5. 
This, however, does not directly allow to check for bifurcations of canonical paths, and, perhaps 
more importantly, requires the recalculation of T at each new u{r]). See |BPS0H iPamOlj for the 
“boundary corrector method” as an approach to avoid the latter, and, moreover, for continuation 
methods in the full t-BVP that for instance also allow the computation of Skiba-curves (cf.^2.1.3|) 
in OD, cf. also iGCF+nSl §7.7-§7.8]. 

As currently p2pOC is based on pde2path, it works most efficiently for spatial 2D problems, while 
ID problems are treated as very narrow quasi ID strips. Presently, pde2path is extended to effi¬ 
ciently treat also ID and 3D problems, based on the package OOPDE www.mathe.tu-freiberg.de/ 
nmo/mitarbeiter/uwe-pruefert/software, Thus, p2pOC will soon provide a genuine ID setting 
as well. 
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